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We study vortex pattern formation in thin ferromagnetic films under the action of strong spin- 
polarized currents. Considering the currents which are polarized along the normal of the film plane, 
we determine the critical current above which the film goes to a saturated state with all magnetic 
moments being perpendicular to the film plane. We show that stable square vortex-antivortex 
superlattices {vortex crystals) appears slightly below the critical current. The melting of the vortex 
crystal occurs with current further decreasing. A mechanism of current-induced periodic vortex- 
antivortex lattice formation is proposed. Micromagnetic simulations confirm our analytical results 
with a high accuracy. 
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I. INTRODUCTION 

The spin-polarized current is a convenient means to 
handle magnetization states of nanomagnets without ap- 
plying of external magnetic fielcP. That is of high ap- 
plied importa nce for constructing purely current con- 
trolled devices^'. One of the effective way to influence 
the film magnetization by the spin-polarized current is 
to use the pillar structure, where the current flows per- 
pendicular to the magnetic filmPHSl Special efforts in this 
way were made to explore the possibility to control the 
properties of magnetic vortexpSHia because the latter is 
a convenient carrier of bit of information. The theoreti- 
cal study in this way is based on the Slonczewski-Berger 
model.EHU 

In this paper we focus on the problem of reg- 
ular pattern formation (vortex-antivortex superlat- 
tice) under the action of strong spin-polarized cur- 
rents, which precedes the saturation. Superstrucures 
of vortices are known from ages of Kelvin's fluid 
vorticesPS Nowadays superlattices of vortices are known 
in superconductivity^, superfluidity^, Bose-Einstein 
condensates (rotating^, nonrotating^, optically dress- 
ing condensate^), and optiaP^H. Vortex-like super- 
structures appear also in magnetism: skyrmion crystals 
were predicted in chiral magnets^, which is now well- 
confirmed experimentally23H2El, vortex-antivortex lattice 
(chirality waves) appears in Kondo lattice modeP^. Re- 
cently we found vortex-antivortex superlattices {vortex 
crystals) in nanomagnets under the action of strong spin- 
polarized current.^ Using micromagnetic simulations we 
found that crystallization precedes a saturation: the 
square superlattices were observed for a range of cur- 
rent densities in immediate vicinity of J c , which is the 
critical current which saturates the magnetization along 
its direction!^ Here we prove theoretically the possibil- 
ity of vortex-antivortex superlattices in ordinary isotropic 
magnetic film. To this end, we build the full theory of 
saturation of a thin ferromagnetic film by transverse spin- 
polarized current. In particular, we show that loss of sta- 
bility of the saturated state leads to appearance of the 



stable square vortex crystals. 

The paper is organised as follows: In Sec. [TT] we de- 
scribe the model and our approach. The linear analysis 
(Sec. Ill) enables us to obtain the value of the saturation 
current J c as function of material parameters and the 
film thickness. The nonlinear analysis (Sec. IV I proves 
the possibility of stable square vortex-antivortex super- 
lattices in pre-saturated regime. All the obtained ana- 
lytical results we check with micromagnetic simulations 
(Sec. [v]). Besides, using the simulations we describe the 
transition from crystal phase into fluid phase which ap- 
pears with the current decrease. 



II. MODEL AND DISCRETE DESCRIPTION 

We consider here a soft magnetic film with thick- 
ness h and lateral size L ^ h. Magnetization of the 
film we model as a three-dimensional cubic lattice of 
magnetic moments M v with lattice spacing a <C h, 
where v = a{v xl V2u v z) with v x ,v y ,v z £ Z is a three- 
dimensional indejP^. In the following we use the nota- 
tions N z = h/a and N xy = L 2 /a 2 for number of mag- 
netic moments along thickness and within the film plane 
respectively. We assume also that the magnetization of 
the film is uniform along thickness. That enable us to 
base our study on the two-dimensional discrete Landau- 
Lifshitz-Slonczewski equationPf^U 



rh n = m n x d£/dm n — jem n x \m T 



(1) 



which describes the magnetization dynamics under in- 
fluence of spin-polarized current which flows perpendicu- 
larly to the magnet plane, along £-axis. It is also as- 
sumed that the current flow and its spin-polarization 
are of the same direction in ([!]). The two dimensional 
index n = a(n x , n y ) with n xi n y G Z numerates the 
normalized magnetic moments m n = M n /\M n \ within 
the film plane. The overdot indicates derivative with 
respect to the rescaled time in units of {4n'-fM s )^ 1 , 7 
is gyromagnetic ratio, M s is the saturation magnetiza- 
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tion, and £ = E/(AttM 2 a 3 J\T z ) is dimensionless mag- 
netic energy. The normalized electrical current den- 
sity j = J I Jo, where J = M 2 \e\h/H with e be- 
ing electron charge and h being Planck constant. The 
spin-transfer torque efficiency function s has the form 
e = ryA 2 /[(A 2 + 1) + (A 2 - l)(m • z)] , where r\ is the 
degree of spin polarization and parameter A ^ 1 de- 
scribes t he mi smatch between spacer and ferromagnet 
resistanc e 1 23 ! 40 !. To simplify representation we omitted 
damping in the equation of motion ([I]), since the role of 
damping is not essential for crystallization of vortices; 
moreover, the saturation current does not depend on the 
damping constant.^ 

The total energy of the system E — E cyi + Ea consists 
of two parts: exchange and dipole-dipole contributions. 
The exchange energy has the form 

Bex = -§ 2 ^ 3imrl ■ m "+«' ( 2 ) 

where n, I are two-dimensional indexes, § is value of spin 
of a ferromagnetic atom, and 3i denotes the exchange 
integral between atoms distanced on I. 
The energy of dipole-dipole interaction is 



E d 



M 2 a b 



-3 



{m v ■ m x ) 



(m„ ■ (A - v)) (m A ■ (A - u)) 



(3) 



where A and v are three dimensional indexes. 
By introducing the complex variable 



Ipr. 



ivrvt 



y/l + m z r 
one can write the Eq. ([I]) in form 



iip r 



d£ 

dip* 



\\4>r, 



l-flV'nl 



(4) 



(5) 



where x = jrj/2 is renormalized current, £ = 1 — A~ 2 and 
ip* denotes the complex conjugation of if). 

It is well known that in the absence of driving (>r = 0) 
the spatially homogeneous state with all moments lying 
in the xy-p\&xie (easy-plane magnetic state) is the most 
energetically favorable state of a thin ferromagnetic film. 
On the other hand, as it is seen from Eqs. |2]), ^ and 
([5]) for large positive k the stationary state of the sys- 
tem corresponds to if) n — or in other words, the system 
goes to the state when all magnetic moments are oriented 
along the z-axis (saturated state). This means that there 
should exist a critical current x c below which the satu- 
rated state loses its stability. Our goal is to study the be- 
havior of the system near threshold of stability of the sat- 
urated state. Near the threshold m z n < 1 and \ip n \ -C 1, 
hence one can expand components of the magnetization 



vector into series in the way similar to the representation 
in terms of the Bose operators! 4 ^ 



rr> 



y _ 4>n - 4>* r 
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mr n = 1 - |</g 2 . 



VP, 
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l^| 2 



0(l^| 5 ) 

o(l^l 5 ) 



(6) 



Substituting §6§ into Q one obtains the equation of 
motion accurate to terms of the third order 



iipr. 



IHtpn ( 1 - 7^\*Pn\ 2 



(7) 



For the future analysis it is convenient to proceed to 
the wave- vector representation using the two-dimensional 
discrete Fourier transform 



(8a) 
(8b) 




with the orthogonality condition 

^2 e i ( fe - fe ')-« — N xy A(k - k') 



(9) 



where k = (k x ,k y ) = ^(l x ,l y ) is two-dimensional dis- 
crete wave vector, l x ,l y € Z, and A(fe) is the Kronecker 
delta. Applying ([8| to the equation Q one obtains equa- 
tion of motion in reciprocal space: 



— llpk = — s h I — — . 



(10) 



where the dimensionless energy of the system can be rep- 
resented as a sum 



£ — £^v £h ~t~ £ 




(11) 



Here the term £° = £° x + £[{ is the harmonic part of the 
energy. It consists of the exchange contribution 



(12a) 



and the dipole-dipole contribution 
-g(kh) 



2 

g(kh) 



1 



(A 



IV^| 2 

x - iky) 2 



k 2 



(12b) 



The nonlinear part of the energy is described by the term 

nl 
d 



£ nl = E 1 ^ + 85 1 , which consists of nonlinear exchange 



cx 

contribution 
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III. HARMONIC APPROXIMATION 



cnl _ I 



E 

fel h2hsk4 



2t(fei,fe 2 )V'fe 1 ^LVife3^fe. 



x A(fei - k 2 + k 3 — fc 4 ) + c.c. 

and dipole-dipole one: 
1 



(13a) 



cnl 



4N r 



E 



x A(fei - fe 2 + fc 3 - fc 4 ) + £(fei)^ fel ^fe 2 Vife 3 Vife 4 (13b) 
X A(fej — fc 2 + &3 + ^4) + c.c. 



The characteristic length 



4irM?a 3 



(14) 



which appears in (|12|) and (13) is so called exchange 



length. We introduced also the following functions 

a(fci,fca) = *i-2(fci-fca), ( 15a ) 

(15b) 

(k x -iky) 2 



B(fci, fca) = g{\k x - k 2 \h) + - 1, 



€(k)=g{kh)- 

, x + e- 
g(x) = 



2k 2 
- 1 



(15c) 
(15d) 



Details of deriving of the Hamiltonian (11)-( |13| i n the 
wave-vector space are placed into the Appendix | A| 

The function £F represents an action of the spin- 
polarized current. It consists of two parts 



1 =5° + 3* 1 , 
with the harmonic contribution 

^° = ^E^ fe 

k 

and the nonlinear part 



(16a) 



(16b) 



Cl,JC 2 ,K3,K4 



4A 2 N X 
x A(fei + k 2 - k 3 - fc 4 ) 



(16c) 



Note that we are interested in a large scale behavior 
of the system and restrict attention to long-wave excita- 
tions. Therefore Eqs. (12), (13) and (15) are written in 
the limit ta « 1. 



First, we discuss solutions of Eqs. (10) in the harmonic 



approximation, since they already capture many essential 
aspects of the problem. By neglecting all nonlinear terms 



in (10), equations for the complex amplitudes tpk and 



ip-k can be written in the form 



-iipk 



2„2 



k 2 l 



g{hk) (k x 



g(hk) 
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iky) 2 



W-h = 



2 k 2 

2n2 , , 9(hk) 
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k Z t - 1 
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iky) 2 



tJ'k 



(17) 



k 2 



(fe)i 



(18) 



The solutions of Eq. (17) have the form form 

Mt) = * +e z+(fc)t , r_ k (t) = *_ e * 

where ^± (fe) are time independent amplitudes and the 
rate constants z±{k) are given by 

z±(k) = -x±x(k), (19) 

where the rate function x{k) is given by 

x(fc) = ^/(l-k 2 e 2 )(k 2 £ 2 +g(hk)-l), (20) 

First of all it should be noted that since x > than 
accordingly to ( [T9| the current plays role of an effective 
damping. That explains why we omitted weak natural 
damping in Eq. (JlJ. That also explains the previous nu- 
merical results, where the saturation and magnetization 
dynamics under the high spin-current influence were in- 
dependent on the damping coefHcientP^. 

Function x(k) is a non-monotonic one which reaches 
its maximum value x r at k = K: 



dk{K) 
dK 







X c 



max 5t(k) 

k 



(21) 



Typical shapes of the rate functions ([20]) are presented 
in Fig. [TJ 

For strong currents when x > x c , we have Re z±(k) < 
for all values of the wave vector k. This means that 
the stationary state of the system is the saturated state 
with m z = 1. However, for x < x c the saturated state 
is linearly unstable with respect to modes tpk with wave 
vectors close to the threshold wave vector K. The corre- 
sponding instability domains for different thicknesses are 
shown in the Fig. [I] as filled regions. 

For each thickness the curve x — x(k), as well as the 
corresponding rescaled curve J(k), separates stable and 
unstable regimes and has a maximum which determines 
the minimal current J c , at which the saturated state re- 
mains stable. So the critical current at which the transi- 
tion to saturation occurs can be determined as 

2M 2 e, 

J c = — i~hx c . (22) 
rjn 
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FIG. 1. Diagram of stability of the uniform state saturated 
transversally by spin current. The regions of instability are 
shown by filling and they are determined by condition *c < x, 
where the the normalized current x is rescaled to the real cur- 
rent density J. Thus for the given current value J < J c one 
has the range [K' , K"\ of the instable wave- vectors. Param- 
eters of material and spin-current were taken the same as for 
simulations (see Sectionjv]). Points show the maximums of de- 
pendencies J(k£) and they determine the saturation current 
for the given thickness. 




h, nm 



FIG. 2. Dependence of the saturation current J c on the film 
thickness. Solid line corresponds to the analytical solution 



obtained from ( 22 I and results of micromagnetic simulations 
(see Section [v]) for different disk diameters D are shown by 
markers. The dashed line demonstrates the parabolic asymp- 
totic for h <C £, see text. 



As one can see from Fig. [T] the saturation current J c 
increases with the increase of thickness. In detail this 
dependence is presented in the Fig. [5] 

Using (20) and (22) one can obtain the following 
h 2 ^ 



asymptotic J c « li 2 ]e\M 2 /{2-qM) for h <C I and J c ~ 
h\e\M 2 j (rjh) for h I though the last one is not achieved 
in the Fig.[2]and it is beyond the limits of applicability of 
the Slonczewski torque in 0. The critical currents ob- 
tained using the micromagnetic simulations appears to 
be in a very good agreement with the theoretical curve. 
Since the present theory is constructed for a film of infi- 



nite lateral size the agreement between simulations and 
the theory is expectedly the best for samples whose thick- 
ness is much smaller than the planar size: h <C D. 



IV. WEAKLY NONLINEAR ANALYSIS 

In this section we prove the stability of structures with 
symmetry C4 which appear in pre-saturation regime. We 
also show that these stable structures are square vortex- 
antivortex superlattices. Thus, our analysis is based on 
the equation ( 10 1 where the nonlinear terms in the Hamil- 



tonian ( 11 ) and in the driving function ( 16 1 are taken into 



account. 

The simplest way to describe the necessary symmetry 
is to restrict ourselves only with four wave-vectors k £ 
{-fCj-, K^, K±, K^} in the wave- vectors space. Here we 
use the following notations 



K 
K 



K(0, 1), 
= K(1,0), 



K 



K(0, 
= K(- 



-1), 
1,0), 



(23) 



where the amplitude K is determined for given thick- 
ness from the linear analysis as following: x c — x(K), 
i.e. K is the wave vector length which maximizes the 
dependence k — x(k). It should be noted that since the 
Hamiltonian £ contains V'fc as well as tp-k then our model 
must contain pairs of vectors (K, —K). It means that 
only structures with even symmetry Cm are possible. We 
focus here on structures with symmetry C4 in order to 
explain the results of the recent numerical experiments^. 

For the future analysis it is convenient to proceed to 
the following notations 



K„ 



(24) 



where a e {t, ->-,4-> •*-}• The value N a in (24) has the 



meaning of the number of magnons with the correspond- 
ing wave vector. Substituting (24) into the equation of 



motion ( 10 ) we obtain the set of eight equations 



N f = 



A 2 U 



(25a) 



xy 



■^/N^N^N^N^ cos($j - $<_>), 



where the other three pairs of equations can be obtained 
by three-time successive rotations of all subscripts by the 
angle 7r/2, and we introduced the notations <&j = $^ + <f>^ 
and $^ = <t>^ + $^ for the sake of simplicity. The 
Hamiltonian in "N — $" -notation being presented as a 
sum of linear and nonlinear parts is the following 



£° 



Mil 



(26a) 



5 



where the linear part reads 



£° = 



i 2 K 2 - 1 



9i 
2 



91 



wVj cos(<i>j) - \/N^ cos($^) 



(26b) 



and the fourth-order nonlinearity has the following form 



l 2 K 2 



1 - 



9i + 92 



£ 2 K 2 + 1 



cos($j - <f\ 



fV2 



91 



9i + 9^2 



(26c) 




cos $ 



r 



iVj- I cos $ 4 



Here we used the analogous notations iVj = Nf + N±, 

= + iV<_, N t = N t N±, = N^N^ and 
<7£ = g(£Kh) to shorten the expressions. 
Using (25a) and (26) one can show that 



df 



(N t - Ni) = -2x(iV t - Ni) 



1 



N t + 2AV 



2A 2 N 



xy 



(27) 



with the corresponding equation for the subscripts ro- 
tated by 7r/2. Taking into account that x > we 
conclude from the Eq. (27) that after period of time 
At = l/(2x) the system achieves a stationary regime 
with N-f = N± and N-y = N^. Consideration of these 



conditions in the stationary form of system (25 ) leads to 
possibility of a solution which satisfy the fol. 
ditions 



lowing con- 



Nf = N± = AT_> = i\r<_ = N, 



(28) 



Under the condition ( 28 1 all four pairs of stationary equa- 



tions of motion ( 25 1 become identical and they obtain the 
following form 



-*(l-H*=-2*(l-§£ 



cos$(l-5,/f)|- =£ 2 K 2 -1 



91 
2 



(29a) 
(29b) 



JIT ( l 2 K 2 - 5 + - 91 + g 2 ) , 



where <yf = NfN xy is density of the magnons, and the 
energy density obtained from the Hamiltonian ( 26 ) reads 



2,yy [2 (£ 2 AT 2 - 1) +.gi(l 



<*)] 



+ 2^ 2 



-£ 2 X 2 + 5 - g 2 - -<?!(!- cos $) 



(30) 




FIG. 3. Magnons density as function of the normalized cur- 
rent for different thicknesses (in units of I): I - 0.5, II - 1, 
III - 2, IV - 4 and A = 2 for all thicknesses. Solid lines show 



the exact numerical solutions of the system ( 29 1 and dashed 



lines correspond to the approximation (31 1. The inset demon 



strates the weakness of influence of the parameter A on the 
exact solution, the data corresponds to the thickness h = 4£. 
Each of the plots is built for the range [x c /2, x c ]. 



Excluding $ from (29) and taking into account that 
jV <C 1 one obtains 



©S + 5x 2 (l 1 



A 2 , 



and then using ( 29b ) one can estimate 



cos $ 



g 

g(Kh) 
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(31a) 



(31b) 



Here we introduced the following thickness dependent 
functions 5 = 2[K 2 £ 2 - 1 + lg{Kh)] and © = [4K 2 £ 2 - 
g(2Kh)]. Magnon density ( |31a ) as well as the corre- 
sponding exact solutions of Eqs. ( p9| are show n in the 
Fig. [3] As one can see, the approximation (31a) is satis- 
factory near the instability threshold. 

Thus we have proved the possibility of a stationary 
structure which is described by the four- wave Ansatz ( 23 ) 
and (24) in the pre-saturated regime. At the same time 



one can see that the parameter A does not influence con- 
siderably the system behavior, see the inset in the Fig. [3] 
The linear stability analysis for the system (25 (shows 



that the stationary solution (31) is stable in the close 
vicinity of the critical current 



< 



< 1, 



see Appendix [B] for details. 

Let us now see how the mentioned structure looks like. 



From (8a I one can obtain the following expression for 



"0-function 



(32) 



G 




L/l 



FIG. 4. The analytically obtained vortex-antivorte x su per- 
lattice. Arrows show distribution of magnetization (331 and 



the corresponding topological density ( 34 1 is shown by gray 
tones. The figure is built for the case h = 41 and x = 0.65x c 
and A = 2 (the required value of $ was determined from ( 29 1 
for the mentioned parameters). 



Varying parameters N a and $ Q one can obtain a wide 
range of different structures from ( 32 1 but under the con- 



ditions ( 28 ) the expression ( 32 1 results exactly the square 



vortex-antivortex superlattice. Indeed, substituting (32) 
into (|6| with taking into account the conditions ( 28 1 we 
obtain in the linear approximation 



2V2T/F 



m y « 2^2^ 



1. 



cos(Kx) sin — + cos(Ky) cos — 

$ $ 

cos(Kx) cos — + cos(Ky) sin — 



(33) 



were the following shift 
performed: x = x + 



of the coordinate origin w; 

* - $<_)/(2.K") and y = y 



($-(■ — <j>j.)/(2iT). Magnetization distribution which cor- 



responds to (|33j) for certain parameters is shown in the 
Fig. [4] by arrows. The topological properties of the sys- 
tem can be characterized by the topological density^ ( or 
scalar chirality densitjE3) v — [d x m x d y m ] ■ m. The 
topological density which corresponds to ( 33 ) reads 



v = -8K 2 yy cos $ sin(Kx) sm(Ky). 



(34) 



The distribution of ( 34 ) is shown in the Fig. [1] by gray 
tones. It resembles chirality waves in Kondo magnets!^ 
Though the model (23 1- (24) results the superlattice 



very similar to one which is observed in the numerical 
experiment, one should point out domain of applicability 
of this model. In p3|-(24) we use only the critical value 



of the wave-vector K instead the whole possible wave- 
vectors in the range [K 1 , K"], where K' and K" bound 



40 



J-JJ, = 0.95 



x« • « • • • • • • • • 



0.6 
0.3 
BO 



-h/l 



FIG. 5. The diagram which determines the range J € 
[J m i n , J c ] for the given geometry sizes (h, L), where the model 
(|23[)-([24| is applicable. In shaded region the model works 
for any currents J < J c . Points correspond to the disks, 
where the superlattices were observed via micromagnetic sim- 
ulations. 



the instability domain for the given current value, see 
Fig. [T] We can restrict ourselves with the single K if 
only size of the system is small enough: 



2tt/L > K" - K'. 



(35) 



Since the domain AK — K" — K' increases with the cur- 
rent decreasing (see Fig. [I]) the condition (35) is equiva- 
lent to J m i n < J < J c , where J min is the minimal current 
at which the condition ( [35] ) . For each of the values of 
thickness h and radius L one can calculate the value of 
Jmin using (20) and (35), the resulting diagram is shown 
in the Fig. [5J Since for the small thicknesses AK -C 1 the 
model ( 23 )— ( 24 ) can be used for wide range of currents. 



V. MICROMAGNETIC SIMULATIONS 

To investigate numerically the process of magnetic film 
saturation under the influence of spin-polarized current 
we used full scale OOMMF 43 micromagnetic simulations. 
All simulations were performed for disk shaped nanopar- 
ticles with material parameters of permalloy: saturation 
magnetization Ms = 8.6 x I0 5 A/m, exchange constant 
A = 13 x 10~ 12 J/m, and the anisotropy was neglected. 
The damping was neglected, because, as it was shown in 
Section [m] the spin-polarized current plays role of an ef- 
fective damping. The mesh cell was chosen to be 3 x 3 x h 
nm. The current parameters r\ = 0.4, and A = 2 were the 
same for all simulations, except some cases mentioned in 
the text bellow. 
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FIG. 6. The superlattice structure obtained using simulations 
in disk with diameter D — 350 nm and thickness h = 20 nm 
under influence of the current J — 32 x 10 12 A/m 2 . Inset 

a) shows the out-of-plane structure of the superlattice, inset 

b) demonstrates in details the magnetization of central part 
of the the disk: arrows correspond to the in-plane magneti- 
zation distribution and out-of-plane component m z is shown 
by color. The superlattice properly is shown in the inset c): 
positions of vortices and antivortices are shown by disks and 
rhombuses respectively. 



In the first stage we obtained the dependence of satura- 
tion current J c on the sample thickness. For this numer- 
ical experiment we chose the nanodisks with three differ- 
ent diameters D — 100, 250 and 450 nm respectivelly, and 
thickness of each of the particles was varied from 0.5 nm 
to 20 nm. As the initial state for a simulation the ground 
state of the particle was chosen: uniform magnetization 
within the sample plane for thin disks (h < 5 nm) and 
vortex state for thicker ones. The spin-current was in- 
creased until the saturation was achieved. As a criterion 
of the saturation we used the relation M z /M s > 0.9999, 
where M z is the total magnetization along the current 
direction. The resulting dependence J c (h) is shown in 
the Fig. [2]by markers. As one can see, for disks with the 
small aspect ratio the micromagnetic simulations confirm 
the analytical results with a high accuracy. The slight de- 
viation from the theoretically predicted curve is observed 
for the case of small disks (see D — 100 nm in the Fig. [2]). 
This is because the presented theory is build for the case 
of an infinite film what corresponds to zero aspect ratio. 

To study the magnetization dynamics in regime J < J c 
we used a disk with diameter D — 350 nm and thick- 
ness h — 20 nm. A spin-current of the certain density 
was sharply applied to this nanodisk, which initially was 
in the vortex ground state. After a few nanoseconds a 
slowly rotating superlattice was formed for case of the 
current value close to the saturation, see the Fig. [6j or 
a fluid-like dynamics of locally ordered vortex-antivortex 
media was observed for cases of lower currents, see Fig. 2 
in the Ref. [38j The Fourier spectrums of the typical 
crystal and fluid structures are compared in the Fig. [7] 
Accordingly to the Fig. [7k) the superlattice is square one. 




-0.4 -0.2 0.0 0.2 0.4 -0.4 -0.2 0.0 0.2 0.4 

*„ 2n/l k,, 2x11 



FIG. 7. Two-dimensional Fourier spectrums of the crystal a), 
and fluid b) structures. Inset a) shows the Fourier transform 
of the function m z (x,y) — (m z ) for the case of the crystal 
structure shown in the Fig. [6] b), where (m z ) is the aver- 
aged m z -component. And the inset b) corresponds to a fluid 
structure obtained for current J = 25 x 10 12 A/m 2 , the other 
parameters are the same as in the Fig. [9] 



To separate crystal and fluid phases and to study their 
properties we performed a series of simulations for a 
range of currents J G [J c /2, J c ] with the current step 
AJ = 0.5 x 10 12 A/m 2 . For a certain value of the cur- 
rent the magnetization dynamics was simulated for 30 ns. 
Starting from the time moment 2 ns we saved the mag- 
netization distribution with the time step 0.2 ns. For 
each of the saved in this way magnetization snap shots 
we found coordinates of all particles (vortices and an- 
tivortices) using the method^ of intersection of isolines 
m x = and m y = 0. To distinguish vortices from an- 
tivortices the winding number of each of the particles was 
calculated as circulation on small circumference centered 
on the particle position. Then for each of the vortices the 
distances to the nearest four antivortices were found (on 
this stage to avoid the boundary influence we consider 
only vortices distanced from the disk center less then a 
half of the disk radius). Then the histogram of the dis- 
tribution of all obtained vortex-antivortex distances was 
build for a certain magnetization snap shot, and finally 
we build the averaged histogram based on all magnetiza- 
tion snap shots for a certain current value. Two examples 
of these averaged histograms are shown in the left column 
of the Fig. [9j The obtained histograms can be well fit- 
ted by the Gaussian f(x) cx exp [—(x — x ) 2 /a 2 ] , where 
f(x) is number of the vortex-antivortex distances which 
are in the interval [x, x + Ax] with Ax = 1 nm being the 
width of the histogram bin. 

For the crystal phase the superlattice constant was 
considered to be a s = 2xq. We found that the superlat- 
tice constant slightly decreases with the current increas- 
ing, this dependence is shown in the Fig. [H] However 
we were not able to determine a s very close to the satu- 
ration current because the components of magnetization 
m x and m y become vanishingly small. The obtained de- 
pendence a s (J) appears to be not very smooth because 
of the stress in the superlattice due to presence of the 
boundary. 
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FIG. 8. Dependence of the superlattice constant as on the 
applied current for permalloy disk with diameter D — 350 nm 
and thickness 20nm. Data were obtained using micromagnetic 
simulations. 



One of the important characteristics, which can be 
extracted from histograms is the cr( J)-dependence. It 
gives a possibility quantitatively separates fluid phase 
and crystal one: in crystals the value of a is small (about 
a few nanometers) and it is weakly dependent on the cur- 
rent J; in the fluid phase the value of a increases fast with 
the current decreasing. To determine the critical current 
Jf c of transition between fluid and crystal phases we fit 
the numerically obtained dependence er( J) by the func- 
tion a = (aJ + b)9(-J + Jf c ) + (aJ fc + b)6(J- J fc ) with 
9{x) being the Heaviside step function and a, b being the 
fitting parameters, see Fig. [9] 



Accordingly to the linear analysis (Section III) the 



parameter A does not influence the saturation current 
J c , and accordingly to the weakly nonlinear analysis of 



the pre-saturated regime (Section IV I the parameter A 
influences very weakly on the dynamics of the vortex- 
antivortex superlattice. Our theory is not able to de- 
scribe the transition between fluid and crystal phases, but 
using the simulations and the methods described above 
we found that the current J/ c , and consequently the cur- 
rent range [J/ c , J c ] of the crystal phase existence, depend 
on the parameter A, see the Fig. [lO] We found out that 
dependence J/ C (A) can be well fitted by the function 



Jfc 







VA2 



1 



(36) 



with p w 7.78 x 10 12 A/m 2 and J° c w 25.92 x 10 12 A/m 2 
being the critical current of the phase transition for case 
A — > oo. For the case A = 1 the superlattice is not 
formed, only the fluid-like dynamics is observed. Also 
the current region of the crystal phase quickly shortens 
when the parameter A is reduced to 1, but for A > 4 the 
the crystal region is approximately constant. 

It should be noted that we do not consider here the 
"gas" phase and the rarefied patterns which were ob- 
served in Ref. [38 at lower currents J< X. 
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FIG. 9. The criterium of separation of fluid and crystal 
phases. In the left column the distributions of the distances 
between the nearest vortices and antivortices are presented, 
solid line shows the Gaussian approximation. The upper 
and lower histograms correspond to the typical fluid-like and 
crystal-like structures respectively. The right plot demon- 
strates dependence of half-width of the mentioned distribu- 
tions on the applied current. All data are obtained from sim- 
ulations for disk with D — 350 nm, h = 20 nm and A = 2. 
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FIG. 10. Phase diagram of the pre-saturated magnetic film 
with thickness h — 20 nm. The transition current J/ c , ob- 
tained from the simulation data (see Section |v| and Fig. |9| is 
shown by points and the corresponding fitting ( |36[ ) is shown 
by the solid line. 



VI. CONCLUSIONS 

We studied theoretically the process of vortex- 
antivortex pattern formation in thin ferromagnetic films 
under the action of strong transversally spin-polarized 
current. We show that there exists a critical (or satura- 
tion) current J c above which the film goes to a saturated 
state with all magnetic moments directed perpendicu- 
larly to the film plane. The critical current strongly de- 
pends on the sample thickness and it is practically inde- 
pendent on the lateral size of the magnet. The saturation 
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current increases with the thickness increasing following 
squared law for thin samples and linear one for thick sam- 
ples. We demonstrate that the stable regular structures 
with symmetry C4 can appear in pre-saturated regime 
and we show that these structures are square vortex- 
antivortex superlattices. Spatial period of the superlat- 
tice slightly decreases with the current increasing. The 
micromagnetic simulations confirm our analytical results 
with a high accuracy. Using the simulations we describe 
the melting of the vortex crystal with the current de- 
crease. 

We show that parameter A which controls the spin- 
transfer torque efficiency, does not modify significantly 
neither the saturation current J c nor the dynamics of 
the vortex-antivortex superlattice. In contrast to this, 
the critical current J/ c which gives the boundary between 
the fluid phase and the crystal one is very sensitive to the 
spin-torque efficiency parameter A: the interval of the 
crystal phase existence [J/ c , J c ] contracts when A — » I 
and it is constant for A> 1. 



Let use perform the Fourier transform ( 8a I with account 
of the orthogonality condition ([9]), which results in 



l k 



5> E 



— ik-2l 1 ^ikal 



x (e- %k21 + e lkil - 2) A(fei - k 2 + fc 3 - fe 4 ) + c.c. 



Now we use the assumption that the magnons whose 
wavelength is of the same order with a are not essen- 
tial for the considered phenomenon, in other words we 
assume that ak <C I. In this case we can expand the 
exponents in exchange energy into series on ak, then 
performing the normalization we finally obtain the ex- 



pressions (12a I and (I3a| for the exchange energy. 



Appendix A: Hamiltonian in the reciprocal space 

Here we calculate the magnetic energy in the wave- 
vector space limiting ourselves by the 4-th order nonlin- 
earity. 



2. Dipole-dipole energy 

In case of the dipole-dipole energy we start from the 
general expression ([3]) . Taking into account that the mag- 
netization is uniform along z-coordinate one can write the 
energy ([3| in formal 



1. Exchange energy 



Let us consider first the exchange energy. Substitut- 
ing the magnetization Q written in terms of %jj into the 
general expression ^ one can write the total exchange 
energy in form £ cx w E®^ + E^, where the linear part 
reads 



E° cx = -S 2 N Z 3, 



n.l^O 



(\ip n \ 2 + \^ n +i\ 2 ) +CC. 



and the corresponding nonlinear part takes the form 



S 2 ^ 



E 



n.l^Q 

|2|„/, |2 



M* n +i(\i> n \ 2 + \^ n +i\ 2 ) 



2|V>n| l^n+ll 1 +C.C. 



Ea 



Mia 



2„6 



E 



A n i (m n mi - 3m* 



+ B nl (m^m? - m v n m\) + C nl (m n m% + m v n m^) 



where the coefficients A, B, and C are the following 



A, 



E 



(K - Vx) 2 + (Xy - Vyf - 2(A Z - V z f 

(A X - V X f - (Xy - Vyf 



2 ^ lA-fl 5 

v*,A 2 1 1 



C n l — 3 2_j 



(Ag - Vx){X V ~ Vy) 
IA-I/I 5 



Vt ,A Z 



Here v — (u x ,p y ,u z ) and A = (X x , X y , X z ) ar e three- 
dimensional indexes while n = (v x ,v v ) and I = (X x ,X y ) 
are the corresponding two-dimensional ones. Substitut- 
ing ^ into the dipolar energy and taking into account 
that A ni = A ln , B nl = Bi n and C n i = C in one ob- 
tains the following expression for dipole-dipole energy 
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E d = El + Ef, where 

M 2 a 6 
2~~ 

71.1 _ 



E° d 



E 



Ef = 



M 2 a 6 



E 



n.l _ 



A„ i (2|^ Il | 2 +^;) (Ala) 

+ Dnil/Jnlpi + C.C. 

a^i 2 (Vi 2 + ;^;) (Aib) 



+ -^D nl \ii n \ 2 ip n il}i + c.c. 



Here the notation D. 



B n i—iC n i was introduced. Now 



we substitute (8a| into (Al) and perform the summation 
over n with taking into account H and finally we obtain 



M 2 n 6 I r 

^=-^El 2i(0) + i(fc) 



+ D{k)\jj k \jj^ k + c.c. 
M?a 6 



(A2a) 



jvinl _ 



- X! I [2i(fei - fe 2 ) + i(fei)] (A2b) 

" fe 1 fe 2 fe3k4V 

X ^fei^fc a ^fes^fc 4 A(fci - fc 2 + fc3 - fcl) 

+ L>(fci)-0 fcl V)fe 2 ^/e3-0fc 4 A(fci - fe 2 + fc 3 + ^4) + C.C 



where functions A(fc) and D(fc) are determined as fol- 
lowing 



A{k) = - E E x ^ + ^ ~ ~ ZXz 



D(k) = 



^ a; 2 - i/ 2 - gjgm 

1 V\ + yf + - z \,)'- 



,5/2' 



(A3a) 

1 1 k 



,5/2 



Form of the functions A(k) and D(k) is not convenient 
for the analysis therefore we perform an approximate 



transition from summation to integration in ( A3 ) . Let 
us start from the function A(k): 

i( fe )=i Dm !&Mh-z) ^ + t: 2 :L ^ xk * +ykv \ 



a* r -*0 

w(r ) 



[x 2 + y 2 + z 2 ] 5 / 2 



where we pricked out the coordinate origin from the do- 
main of integration w(ro), see the Fig. 11 and we also 
used the relation 

h h h 

j dz J dz'F{\z - z'\) = 2 J(h- z)F(z)dz. (A5) 
00 



I 




w(r ) 


h \ 











FIG. 11. Cross-section of the film. The filling shows domain 
of integration w(ro) used in the calculation of function A(k) 
in JA4l. 



Separating the region of integration w(ro) into parts I 
and II (see Fig. 11) and performing the change of vari- 
ables (a;, y) — p(cos\, sin^) we can represent the func- 
tion A as a sum A(k) = {Aj + An) /a 4 , where 



Ai 



2iT =lun Q Jdz(h-z) Jdp P [ P + * Z ]5/2 Mpk), 
Po 



(A6a) 

po 00 

2? = P ^oJ dz{h - Z) J d ^ [p 2 + ~j] 5/2 J oW 



(A6b) 

Here we performed integration over x € [0; 27r] us- 
ing the relation e^O* cos x+ fcH sm X )d x = 2tt J (fcp), 
where Jo (a;) denotes zero-order Bessel function of the 
first kind. Direct integration in (A6a) with the conse- 



quent limit calculation results A//27T — —hg(kh), where 
g(x) = (e x + x — l)/x. Change of variables p — > pop and 
z p {) z allows us to get rid of the po in the integration 
limits: 

~ = Urn f dz(h-zp ) [ dpp 9 2 * J Q (pp k) 
2tt po->-oj J [p z + z- ! \ t >> z 

(A7) 

Calculation the limit in (A7| with the subsequent inte- 

(A8) 



(A3b) gration results Ajj/2n = h2/3, so finally 



A(k) 



2-Kh 

-.4 



a' 



g{kh) 



Performing the same transition to the polar coordi- 
nates (p, x) we represent (A3b) as following 



D(k) 



1 



dz(h — z) 



p 3 dp 



2tt 



[p 2 + Z 2 fl 2 





dxe 



-i{pk+2 X ) 



(A9) 

where pk = p(k x cosx + k v sinx). The direct integration 

2tt 



jA9| ) using the relation J d X e^ x cos x+nx ~> = 2ni n J n (x) 
results 



{) = -^r 9 ^ — k 2 • 



(A10) 
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Substituting now ( A8 ) and ( A10 1 into ( A2 1 we obtain the 



expressions ( 12b I and ( 13b ) for the dipole-dipole energy. 



Appendix B: Stability of the vortex-antivortex 
lattice solution 



Here we consider stability of the stationary solu- 
tion (28), (31) of the system (25). As it was sown in 



the main text after period of time r = 1/ (2k) the so- 
lutions of Eqs. (25) satisfy the conditions N-f — N± and 
N_> = . Using these conditions and introducing vari- 
ables Ni = JVf = Ni, N 2 = = N^, $i = $ t + $4. 
and $2 = + one can reduce the system of eight 



equations ( 25 1 to the system of four equations 
9£ 



Ni = -w t - F * ■ 

1 dN t 1 ' 



(Bl) 



i = 1, 2. 



Here the forces of the spin-current acting has the form 

3A, + ANv t 



Ff =2x|iV, 
2k 



1 - 



2A 2 N 



A 2 K 



cos($i - $ ? ) 



A 2 K. 



iV-sin^-^), 



(B2) 



where the notation 1 = 2 and 2 = 1 is used. Hamiltonian 



( 26 1 in terms of the new variables takes the form 
£ = £° + £ nl . 



where the linear part ( 26b ) reads 



£° = 2(N X + N 2 ) (kH 2 + f - l) ~ 
- gi (Ni cos $i - N 2 cos $ 2 ) 



(B3a) 



(B3b) 



and the corresponding nonlinear part ( 26c 
gal = 



is 



xy U=l,2 



■K'r-n9i-9») 



(B3c) 



+ 2A^(cos($, ; - + K 2 £ 2 - ^ - g^)+ 

+ 2 - gi - 5v3 )] + ~ gi (N 2 cos - N 2 cos $ 2 ) 

giNiN 2 (cos&i - cos$ 2 ) >■ 



Now we linearize the system (Bl) against a stationary 
solution v° ={N?, N°, $°, $°}~ 



Mv. 



where v is small deviation from the solution v° and 4x4 
matrix M can be presented in the following block form 

/ M NN M iV*\ 
M — I ]yj*JV jy[** I ) (Boa) 

\ / v— v° 

where the components are the following 2x2 matrixes 



M NN 



M 



d 2 £ 


OF? 


dNidNj 


dN 3 


d 2 l 


OF? 


dNid<5> 3 


d$j 


d 2 l 


dFf 



(B5b) 



d$idNi dNi 



M 



a 2 £ 



9^ 



1,2. 



Necessary and sufficient condition for stability of the 
solution v° is negativity of real parts of all eigenval- 
ues Ai, A 2 , A3, A4 of the matrix M. After the straight- 
forward calculation of ( B5b[ ) we substitute the solution 
Ni = N 2 = N, $1 = $ and $ 2 = $-7r what corresponds 



to the conditions (28). Then after the straightforward 
calculation of the eigenvalues of M we exclude $ using 
( 29 1 and finally we obtain 



Ai 



A, = 



y/2. 



y/2. 



s]2*l - >r 2 - 
^2kI ~k^ + 



2Jf 
x c 

2Jf 
x c 

2Jf 

1JV 

K r 



3® 
5® 
3® 
[3®- 



3k 2 1 



hK 1 



5x; 



1 

A 2 
1 

" A 2 
1 

" A 2 



(B6) 

where only the terms linear with respect to ,jV are saved. 
It should be noted that for JV = (what corresponds to 
the saturated state) the condition Ai < in ( B6 1 is equiv- 



alent to the condition of s tabilit y of the saturated state 
k > k c . Substituting (31a I into (B6 1 and considering the 



indefinitely small deviation from the saturation current 
k = k c +5 one obtains the following linear approximation 
of ( B6 1 with respect to deviation S 



Ai 



A 



26 



3® 



1 

A 2 



45 
28, 
A 4 = -2 



3© 



3® + bx 2 c (1 
3k 2 (1 - i) 



A 2 - 



(B7) 



A 



26 



3® + 5x 2 



3® + 5>r 2 (1 



A 2 - 



Taking into account that 3® > and A > 1 one can 
conclude that all Ai < only for 8 < 0. That means 
that the obtained lattice solution (28), (pll) is stable for 



(B4) infinitely small decrease in current from the critical value. 
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